library(tidyverse)
library(redist)
library(dataverse)

# Functions from fifty-states
source("R/00_setup.R")
source("https://raw.githubusercontent.com/alarm-redist/fifty-states/main/R/summary_stats.R")
source("https://raw.githubusercontent.com/alarm-redist/fifty-states/main/R/validate.R")

# Data ------
# raw dataverse file of map object
map_2020_dv <- get_dataframe_by_name(
  filename = "OH_cd_2020_map.rds",
  dataset = "10.7910/DVN/SLCD3E",
  server = "dataverse.harvard.edu",
  .f = read_rds
)

map_2020 <- map_2020_dv

# Setup -----
set.seed(2020)

N <- 30000 # simulations
# removed custom columbus constraint that was here

# Simulate ------
plans <- redist_smc(map_2020, N,
  runs = 2,
  counties = county,
  pop_temper = 0.04,
  seq_alpha = 0.95, verbose = TRUE
) |>
  pullback(map_2020) |>
  group_by(chain) |>
  filter(as.integer(draw) < min(as.integer(draw)) + 2500) |> # thin samples
  ungroup()

plans <- match_numbers(plans, map_2020$cd_2020)


# Compute summary statistics -----
splits_mat <- apply(as.matrix(plans), 2, \(x) tapply(x, map_2020$county, n_distinct)) - 1L
nd <- attr(map_2020, "ndists")

plans <- add_summary_stats(plans, map_2020) |>
  mutate(
    splits_1 = rep(as.integer(colSums(splits_mat == 1)), each = nd),
    splits_2 = rep(as.integer(colSums(splits_mat == 2)), each = nd),
    splits_3 = rep(as.integer(colSums(splits_mat >= 3)), each = nd)
  )

write_rds(plans, "data/ohio_remove-constraints.rds")


# plot results------
summary(plans)
validate_analysis(plans, map_2020)
